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Thermal transport is an important energy transfer process in nature. Phonon is the major energy carrier for 
heat in semiconductor and dielectric materials. In analogy to Ohm's law for electrical conductivity, Fourier's law 
is a fundamental rule of heat transfer in solids. It states that the thermal conductivity is independent of sample 
scale and geometry. Although Fourier's law has received great success in describing macroscopic thermal 
transport in the past two hundreds years, its validity in low dimensional systems is still an open question. Here 
we give a brief review of the recent developments in experimental, theoretical and numerical studies of heat 
transport in low dimensional systems, include lattice models, nanowires, nanotubes and graphenes. We will 
demonstrate that the phonon transports in low dimensional systems super-diffusively, which leads to a size 
dependent thermal conductivity. In other words, Fourier's law is breakdown in low dimensional structures. 



INTRODUCTION 

The conduction of heat is one fundamental energy trans- 
port mechanisms in nature. In addition to electron and pho- 
ton, phonon- the heat pulse through lattice, also carries and 
processes energy, which has broad applications for heat con- 
trol/management in the real world. Actually, for non-metallic 
materials, phonon is the most dominated heat energy carrier. 
Its contribution to heat conduction is much larger than those 
from electrons and photons. Traditionally, the phenomenon 
of thermal transport is believed to follow the Fourier's law of 
heat conduction 



J = -kVT, 



(1) 



where J is the heat flux in the system, VT is the gradient of 
temperature. The conductivity k is a geometry-independent 
coefficient which mainly depends on the composition and 
structure of the material and the temperature. This empiri- 
cal law has received great success in describing macroscopic 
thermal transport and is widely accepted as a general truth. 

In the past decades, low dimensional nano scale systems 
have been extensively studied due to their promising po- 
tential applications for future electronic, optoelectronic, and 
phononic/thermal devices. In addition to the electrical and op- 
tical properties, the thermal properties of nanoscale materials 
are also important but less studied compared with electronic 
and optical properties. The thermal properties in nanostruc- 
tures differ significantly from their bulk counterpart because 
the phonon characteristic lengths are comparable to the char- 
acteristic length of the nanostructures. However, it is still an 
open and much debated question whether the Fourier's law is 
valid in low dimensional systems. A rigorous proof for this 
empirical law from microscopic Hamiltonian dynamics is still 
absent even though it is already two-hundred-year old. Here 



"invalid" means the thermal conductivity k will depend on the 
system size. 

Such size dependent thermal conductivity has been ob- 
served in theoretical models, such as the harmonic chains 
[1, 2], the FPU-/? model [3, 4] and the hard point gas model 
[5-7]. The discovery of this anomalous behavior in general 
low dimensional models has then inspired enormous research 
studies. It is found that in these models, the thermal conduc- 
tivity typically diverges with system size: 



for ID models and 



k ~ log L 



(2) 



(3) 



for 2D models in thermodynamic limit, which predicts ex- 
traordinarily high thermal conductivity for low dimensional 
systems with finite but large sizes. Therefore, not only it 
is a fundamental demand for the development of statistical 
physics to understand normal and anomalous heat conduction 
in low dimensional systems, but it is also of great interest from 
the application point of view, since the achievement of modern 
nano fabrication technology allows one to access and utilize 
one dimensional (ID) and two dimensional (2D) structures 
with sizes in the range of few nanometers up to few hundred 
nanometers. 

In recent years, much effort has been devoted to study the 
anomalous heat conduction in low dimensional nanostructures 
by numerical simulations such as [8] and [9]. These studies 
provide evidence that low dimensional nanostructures are very 
promising platforms to verify fundamental thermal transport 
theories. Inspired by these theoretical and numerical stud- 
ies, it is exciting that length dependent thermal conductiv- 
ity has been observed in carbon nanotubes and boron-nitride 
nanotubes [10]. Unfortunately, an experimental investiga- 
tion of the quantitative size dependent thermal conductivity in 
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nanoscale is exceptionally challenging. The difficulty lies pri- 
marily with the technique associated with measuring temper- 
ature distribution and heat flux of low dimensional systems, 
with the added complexity of accurately varying the size of the 
investigated object. Due to these challenges in experimental 
measurement, the experimental data need to be complemented 
by theoretical study on a microscopic footing. 

Generally, heat transport can be studied in the aspects of 
two different yet deeply related phenomena. One is to study 
the heat flux in a non-equilibrium steady state to get the heat 
conductivity as we normally do in simulation or in experi- 
ments, in which temperature biases are applied to the system. 
The other approach is to study the heat diffusion in a non- 
equilibrium transient process in which the system is firstly ex- 
cited away from equilibrium and then relax to equilibrium. 
Both approach should reveal some properties of heat conduc- 
tion of the system. Moreover, the connection between thermal 
conductivity and heat diffusion, especially for the anomalous 
cases, is still an open question and attract wide research inter- 
ests. 

As heat conduction in low dimensional systems is a fast 
growing area, combination of experimental, theoretical and 
numerical investigations are indispensable to speed up its de- 
velopment. In this articles, we would like to give a review 
on the recent development of heat transport in low dimen- 
sional systems, from both experimental and theoretical point 
of view. The rest of the article is organized as follows: section 
introduces the experiments on quasi- ID nanostructures and 
2D graphene that display the anomalous size dependent ther- 
mal conductivity. Section discusses the numerical studies on 
thermal conduction of nanotubes, nanowires and graphenes, 
which would shed light on understanding the phenomenon of 
size dependent thermal conductivity from a theoretical point 
of view. Section is devoted to the general theories on anoma- 
lous heat conduction and anomalous energy diffusion in low 
dimensional systems. Finally, in section and , we present a 
few conclusions and a short outlook. 

There exists a large amount of literatures and studies on dif- 
ferent aspects of thermal property in low dimensional systems. 
For a comprehensive review on thermal property of nano ma- 
terials, please refer to these articles Refs. [11-13]. There are 
also a few reviews talk about anomalous heat transport in low 
dimensional systems from the viewpoint of fundamental sta- 
tistical physics, such as Ref. [14] and [15]. Due to the limit of 
space, we only address the most fundamental physical issue, 
i.e., the violation of Fourier's law in this article. 

EXPERIMENTAL OBSERVATION OF ANOMALOUS HEAT 
CONDUCTION IN REAL NANO MATERIALS 

Anomalous Heat Conduction in Quasi ID Nanostructures 

When the dimensions of materials shrink to the nanoscale, 
thermal transport properties can be very different from those 
of their counterpart bulk materials. Recent experimental in- 



vestigations in (quasi) ID structures have revealed two oppo- 
site intriguing properties: unusual high thermal conductivities 
and significantly suppressed thermal conductivities. The un- 
usual high thermal conductivities are attributed to the unique 
crystalline structures of the ID materials, such as carbon nan- 
otubes [16], and stretched polymer nanofibres [17], whereas 
the significantly suppressed thermal conductivities are due to 
an increased phonon-boundary scattering, which has been ob- 
served in various nanowires, such as Si [18-20], Si/SiGe su- 
perlattice [21], ZnO [22], Bi [23], etc. More recently, it has 
been observed that the thermal conductivity can be further 
suppressed by the coherent phonon resonance effect beyond 
the phonon-boundary scattering limit in Ge-Si core-shell ul- 
trathin nanowires [24]. 

Although significant experimental efforts have been de- 
voted to study the diameter-dependence of the thermal con- 
ductivity to understand the phonon-boundary scattering mech- 
anisms in the nanostructures, the study of length dependence 
of the thermal conductivity to address the validity of Fourier's 
Law in ID nanostmcture remains rare. The primary techni- 
cal challenges lie in the problem of suspicious thermal con- 
tact resistance and the difficulty in guaranteeing the diame- 
ter exactly the same when varying the length of the samples 
in different measurements. In 2008, Chang et al. developed 
a sequential multiprobe method that can establish the devia- 
tion from Fourier's Law behavior in one-dimensional nanos- 
tructures, without suffering from the problem of suspicious 
contact resistance [10]. A suspended microelectromechani- 
cal system (MEMS) device was used to measure thermal con- 
ductivity of carbon nanotubes (CNTs) and boron-nitride nan- 
otubes (BNNTs) as a function of length. Unlike the traditional 
method using a movable probe as a local probe, Chang et al. 
deposited a series of thermal contacts to vary the length along 
the same sample, thus ensuring the diameter of the sample 
unchanged. 

Fig. 1 (a) shows experimental procedure, and the inset in Fig. 
1(a) shows the nanotube bonded to the test fixture by Pt-C 
composite using electron-beam-induced deposition. The right 
contact attaches the nanotube to the top of a preformed verti- 
cal "rib" on the right thermal pad; this elevation ensures that 
the nanotube is fully suspended between the two contacts. In 
order to vary the length of the suspended segment, a series of 
additional thermal contacts, labeled as 1 to 5 in Fig. 1(a), was 
deposited inward of the original right-hand contact using Pt-C 
composite. 

The thermal conductance of the nanotube as a function of 
length was measured by the method previously described by 
L. Shi et al. [25]. In brief, the Pt loop on one of the suspended 
thermal pads acts as a heater , while the other loop acts as a 
sensor (Fig. ??a). Heat loss from the heater is through the 
low-thermal-conductance leads suspending the thermal pad, 
and through the sample bridging to the sensor. By measur- 
ing the temperature rises at the heater and the sensor, taking 
into account the heat loss through its suspending leads, the 
thermal conductance of the sample can be calculated. From 
the analysis of nanotube thermal conductance versus sam- 
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The different length behavior might be due to the difference 
in isotropic disorder between BNNTs and CNTs. It is worth 
pointing out that the observed f3 on multi-walled CNTs dif- 
fers from the theoretical predictions on single-walled CNTs 
[8]. The possible physical origin for this discrepancy might 
be the inter-shell phonon scattering in multi-walled CNTs. 

In addition to the experimental observation of the viola- 
tion of the Fourier's Law in both BNNTs and CNTs, the 
anomalous heat conduction was also experimentally observed 
in polymer nanofibres by Shen et al. [17]. Bulk polymers 
are generally regarded as thermal insulators because they have 
very low thermal conductivity on the order of 0. 1 W/mK. In- 
spired by the theoretical work which suggests that individual 
polyethylene chains can have extremely high thermal conduc- 
tivity [26], Shen et al. fabricated high-quality ultra-drawn 
polyethylene nanofibres with diameters of 50-500 nm and 
lengths up to tens of millimeters. Shen et al. used a sensi- 
tive bi-material AFM cantilever in combination with a micro 
thermocouple to measure the thermal conductivities of the in- 
dividual nanofibres, as shown in Fig. 2a. It was found that 
the nanofibres have an extraordinary high thermal conductiv- 
ity which increases with increasing draw ratios (Fig. 2b). The 
highest measured thermal conductivity (104 W/mK) is about 
three times higher than that of micrometer-sized fibres and 
300 times that of bulk polyethylene ( 0.35 W/mK). A value 
of 104 W/mK is higher than many metals, including Pt, Fe 
and Ni. The high thermal conductivity was attributed to the 
restructuring of the polymer chains by stretching, which im- 
proves the fibre quality toward an ideal single crystalline fibre. 
Very recently, the similar anomalous high thermal conductiv- 
ity is also observed in spider silk [27]. 



FIG. 1. (a) upper: SEM image of a thermal conductivity test fixture 
with a BNNT after five sequences of Pt-C deposition. The numbers 
denote the nth deposition. The inset shows the SEM image after the 
first Pt-C deposition. The arrow denotes the preformed 'rib' for sus- 
pending the BNNT. Lower: Two circuit models for analyzing the data 
taking into account of the contact resistance. Rs(n) and Rc(n)denote 
the sample resistance and the contact resistance at each deposition, 
respectively. (b) Normalized thermal resistance vs normalized sam- 
ple length for CNT sample (solid black circles), best fit assuming 
/3=0.6 (open blue stars), and best fit assuming Fourier's law (open 
red circles), (c) Normalized thermal resistance vs normalized sam- 
ple length for BNNT sample (solid black diamonds), best fit assum- 
ing /3=0.4 (open blue stars), and best fit assuming Fourier's law (open 
red circles). For further details see in [10]. 



pie length, taking account into the finite contact resistance 
between the nanotube and the thermal pads (Fig. 1(a)), it is 
clear that thermal conductivity of nanotube does not follow 
the Fourier's law, regardless of whether L ^> A (A is phonon 
mean free path), as shown in Fig. 1(b) and Fig. 1(c). Chang 
et al. . found that the thermal conductivity diverges with tube 
length as k a L@(L is the length of the nanotube; (3=0 corre- 
sponds to Fourier's Law). The value of (3 ranges from 0.6 to 
0.8 for CNTs, whereas f3 ranges from 0.4 to 0.6 for BNNTs. 



Anomalous Heat Conduction in Graphene 

Although significant progress has been made for ID sys- 
tems, the study of heat conduction in 2D systems is still in 
its infancy due to the lack of true 2D materials, lack of 2D 
phonon transport theory and short of computational power . 
The discovery of graphene has changed this scenario [28]. 
The first experiment, based on Raman spectroscope measure- 
ment (Fig. 3), was carried out by UC riverside group in 2008 
[29-31]. This experiment was performed on suspended single 
layer graphene (SLG) exfoliated from HOPG. The room tem- 
perature thermal conductivity k was measured to reach a high 
value of —4800-5300 W/mK, exceeds that of bulk graphite 
and carbon nanotube [16]. Such a high thermal conductiv- 
ity is related to the long phonon mean free path of graphene, 
which is calculated to be around 750 nm at room temperature 
[29]. 

The thermal conductivity in fully supported graphene is 
however much smaller. The measurement of Si02-supported 
graphene revealed a value of —600 W/mK at room temper- 
ature, which is due to backside scatterings and the flexural 
phonons leakage into the substrate [32]. By introducing an 
exact numerical solution for Boltzman's transport equation, 
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FIG. 2. (a) Schematic of experimental set-up used to measure the 
thermal properties of a single ultra-drawn nanofibre. The thermal 
sensor is a silicon nitride AFM cantilever coated with a 70-nm gold 
film. A laser beam (wavelength, 650 nm; output power, 3 mW) is 
focused on the tip of the cantilever and reflected onto a bi-cell pho- 
todiode. The nanofibre drawn from the AFM cantilever is loosely 
suspended between a micro thermocouple and the AFM cantilever, 
(b) Thermal conductivities of three samples vs their corresponding 
draw ratios. For further details see in Ref. [17]. 



Seol et al. claimed that the thermal conductivity of suspended 
graphene should be around 5 times larger, which is ~ 3000 
W/mK at room temperature and is consistent with the results 
from UC riverside group. 

Several following Raman based measurements, performed 
by other independent groups, confirmed the high thermal con- 
ductivity both in exfoliated and CVD graphene with corbino 
geometry [33-36]. Cai et al. reported that the thermal con- 
ductivity of suspended CVD graphene exceeds ~2500W/mK 
near 350K and decreases to around ~1400 W/mK at about 
500K [33]. Lee et al. found that the value in suspended exfo- 
liated graphene is ~ 1800 W/mK near 325K, and decreases to 
~700W/mK at around 500K [35] . Another group repeated the 
similar experiment on exfoliated graphene and found a value 
of ~600W/mK [36]. The differences can be explained by the 
difference in the actual temperature in graphene under laser 
heating, optical absorption, laser spot size, and the well-know 
uncertainty in the Raman based measurements [37]. Despite 
the variation in the measured values, the thermal conductivity 
in both suspended and supported graphene is believed to be 
much larger than that in silicon and most of the metals. 

Apart from Raman based measurement, typical thermal- 




FIG. 3. Schematic of the Raman based thermal conductivity mea- 
surements. The focused laser light creates a local hot spot and the 
temperature rise is measured by Raman spectroscope. For further 
details see Ref . [29]. 




FIG. 4. (a)-(b) SEM images of the thermal bridge for suspended and 
supported graphene samples; scale bar: 5 /im. (c)Zoom in image of 
the thermal bridge, corresponding to (a); scale bar: 1 fim. (d) Equiv- 
alent thermal circuit of the device, KsiNx is the thermal resistance 
of the nitride platform for supported samples. For further details see 
Ref. [38], 



bridge method, using the MEMS membranes mentioned 
above, was also used to measure the thermal conductance of 
graphene by the group from National University of Singapore 
[38, 39]. Fig. 4a and 4b show scanning electron microscope 
(SEM) images of the suspended thermal-bridge used for the 
measurements of suspended and supported graphene, respec- 
tively. For details of the fabrication process and measurement 
accuracy, please refer to L. Shi et al. [25]. 

Suspended CVD SLG was measured with length of 300 
nm by Xu et al. using the thermal-bridge method [39]. o7A 
(thermal conductance per unit cross section area) was found 
to reach a high value of around 1.8 x 10 5 T 1 53 W/m 2 K below 
100K (Fig. 5). Mingo and Broido [40] have predicted that the 
thermal conductance will approach the ballistic limit in clean 
devices, with calculated a/A to reach 6xl0 5 T 15 W/m 2 K. 
The experimentally observed values are within 30% of the 
predicted ballistic thermal conductance in graphene, demon- 
strate that the thermal transport in suspended micron graphene 
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FIG. 5. Thermal conductance per unit cross section area a/ A as a 
function of temperature. The measured data are approaching the ex- 
pected ballistic limit (black dashed curve). For further details see 
Ref. [39], 



is indeed in the ballistic regime. Furthermore, the temperature 
dependence of k can be well fitted with k =bT" (a ~ 1.5). 
This observed temperature dependence of k(T) and the high 
value in a/ A clearly identify the dominant phonons contribu- 
tion to the thermal conduction in clean graphene. Mingo et 
al. have argued that the ZA phonons carry most of the heat in 
SLG; at low temperatures, the out-of-plane ZA modes are pre- 
dicted to lead to a T 15 behavior of the thermal conductivity 
[41]. 

In another experiment, the room temperature thermal con- 
ductivity was measured to be around lOOOW/mK in sus- 
pended SLG with 5 /im in length and 1.6 /jm in width [42], 
which is 40% to 50% larger than that obtained using the 
same measurement technique in suspended exfoliated bilayer 
graphene of similar sample size [43]. This is expected due to 
the lack of interlayer coupling in SLG [31]. While the Raman 
based measurements in C VD/exfoliated SLG show an average 
value which is much larger [29, 33, 34], a direct comparison 
is challenging raised from the uncertainty of Raman measure- 
ments [37], and, more likely, the differences in sample geom- 
etry. 

Efforts have been made to study the size effect of thermal 
conductivity in SLG. Chen et al. studied the thermal conduc- 
tivity in suspended CVD SLG of corbino geometry with dif- 
ferent diameters ranging from 2.9 /im to 9.7 /im using Raman 
method [34]. The measured thermal conductivity randomly 
varies from ~ 2600 W/mK to ~ 3100 W/mK with diameter 
changes, showing no sample sized dependence. The authors 
attributed this to a relatively large measurement uncertainty as 
well as grain boundaries, defect or possible polymer residues. 

Alternatively, it was found experimentally that the ther- 
mal conductivity of multi-layer graphene changes with sam- 
ple size [38]. In this experiment, both the suspended and 
supported multi-layer graphene were measured by using the 
suspended thermal-bridge configuration mentioned above. A 
room temperature thermal conductivity of 1250 W/mK was 
obtained for a 5 /itn-long supported flake (Fig. 6). Experimen- 
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FIG. 6. Thermal conductivity of multi-layer graphene vs. tempera- 
ture. The lengths of samples SI, S2 and S3 are 5 /im (three layers), 
2 /im (five layers) and 1 /im (three layers), respectively. All the sam- 
ples have the same width of 5 /im. Sample S 1 and S3 are supported 
by a thin SiN^ membrane and S2 is suspended. For further details 
see Ref. [38]. 



tal results show that thermal conductivity depends strongly on 
the graphene size and is lowered by 85% when the length 
of the flake reduces from 5 /im to 1 /mi (SI and S3 in the 
same Figure). This length dependence is consistent with the 
calculation in SLG that thermal conductivity is predicted to 
be length dependent even when graphene flake is as long as 
100 /tm, due to the emergence of the long wave phonons 
when sample size increases [44]. More interestingly, theories 
had envisioned the length dependent behavior in 2D systems 
before the thermal conductivity in graphene was carried out 
[15, 45]. A detailed discussion on the size dependence of the 
thermal conductivity in a 2D Fermi-Pasta-Ulam lattice was 
carried out and a ~ logL behavior is proposed in a wide range 
of L and W/L ratio [46]. 



ATOMISTIC SIMULATIONS OF HEAT CONDUCTION IN 
NANO MATERIALS 

It is obvious that the systematic applications of nanomate- 
rials based devices will be greatly accelerated by a detailed 
understanding of their material property. For nanoscale ma- 
terials, the fundamental question one may ask is whether 
Fourier's law is still valid. The question is not trivial, since 
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nanoscale structures is of finite number of atoms and far 
from the thermodynamic limit. Due to the challenge in 
nanoscale experimental measurement of temperature and heat 
flux, atomistic simulations have an important contribution to 
the development of this area. In this section, we focus on 
the simulation investigations about the thermal conduction 
and diffusion in low dimensional nano materials. Because 
a large variety of studies on thermal transport of nanoscale 
materials have been done in the past decade, here we only ad- 
dressed the most fundamental aspects of validity of Fourier's 
law in nano materials, in particular, in nanotubes, nanowires 
and graphene. 



Anomalous Heat Conductivity and Energy Diffusion In 
Nanotubes and Nanowires 



Carbon nanotube (CNT) is one of the promising nanoscale 
materials come to the spotlight of research after it was discov- 
ered in the 1990s [47]. In addition to electronic and optical 
properties, thermal property of CNTs has attracted more and 
more interests. It was found experimentally that at room tem- 
perature the thermal conductivity of a single CNT is about 
3000 W/mK [16]. Recently, thermal contact resistance of 
carbon nanotube junctions have been explored experimentally 
[48]. Both electron and phonon can be heat carrier. Yamamoto 
et al. has demonstrated that even for the metallic nanotubes, 
the electrons give limited contribution to thermal conductiv- 
ity of single-walled CNTs (SWNTs) at low temperature, and 
with increase of temperature this part decreases quickly [49]. 
It is found that the phonon mean free path of CNT can ex- 
ceed the characteristic length of the structural ripples. This 
makes CNT an ideal phonon waveguide which can have supe- 
rior phonon transport property to its counterpart of photon in 
optical fibers [50]. Thus SWNT is ideally suited for molecular 
dynamics (MD) investigations of the phonon thermal conduc- 
tion law in low dimensional system. By using MD simula- 
tions, it was found that the thermal conductivity of SWNT 
diverges with the length of system as k ~ L 13 , where the ex- 
ponent /3 depends on the temperature and SWNT diameters, 
and the value of j3 is between 0.12 and 0.4 [51] and between 
0.11 and 0.32 [52]. As discussed in above section, this length 
dependent thermal conductivity has been verified experimen- 
tally, although the absolute value of (3 varies. 

To understand the physical mechanism for the length de- 
pendent thermal conductivity in SWNT observed both exper- 
imentally and theoretically, Zhang and Li investigated the vi- 
bration energy diffusion along SWNT [8]. In the calculation 
the Tersoff empirical potential is used to derive the force form, 
which has been widely used in the study of heat conduction 
in nanostructured carbon systems. They first thermalize the 
CNT to a temperature T, then a packet of energy (heat pulse) 
is excited in the middle of the tube to study how it spreads 
along the system. To suppress statistical fluctuations, an av- 
erage over 1000 realizations is performed. Fig. 7(a) shows 
the representatives of pulse propagation along the lattice at 
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FIG. 7. (color online). Energy diffusion in carbon nanotubes. (a) 
Representatives of pulse propagation along the SWNT at 2K. (b) The 
snapshots of heat pulse diffusion in SWNT at 300K. (c) and (d) En- 
ergy diffusion in SWNT at 2K and 300K, respectively. For further 
details see in Ref. [8]. 



2K. The pulse diffusion profiles at room temperature are very 
similar to those at low temperature. In the SWNT, one sin- 
gle peak expands with time, and it is reflected back from the 
two ends after the wave front reaches the tube ends, no peak 
splits are observed. Indeed, there exist three time periods in 
the diffusion process. At the start period, the pulse spreads 
between the two nearby layers, which is denoted as "initial 
period." Then, the pulse diffuses along the tube, this is "dif- 
fusion period"; finally, the wave front hits the boundary and 
is reflected back and transport in opposite direction, which 
is "reflection period." The snapshots of representative pulse 
propagation along the CNT at room temperature is shown in 
Fig. 7(b). It is clear that in CNT, one peak expands with time, 
this pulse spread can be described quantitatively as below: 



£Jflj(*)-ifr]M*)-ri(0)) 2 



(4) 



here Ei(t) is the energy distribution of atom % at time t, r,(0) 
is the position of energy pulse at t = 0. 
The averaged energy profile spreads as: 



a 2 {t) = 2Dt a , with < a < 2, 



(5) 



where (.) denotes the ensemble average over different realiza- 
tions. 

In Fig. 7(c) and (d) we show (cr 2 (t)) versus time in double 
logarithmic scale, so the slope of the curve gives the value of 
a. It is clearly seen that the slope is two at 2K. This means 
that phonon energy transports ballistically at low tempera- 
ture. This can be understood from the Taylor expansion of 
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the atomic interaction potential by keeping up to the second- 
order term. At low temperature, the vibrations of atoms are 
very small, the Tersoff potential can be approximated by a har- 
monic one, thus corresponding to a ballistic transport. How- 
ever for SWNT at room temperature, energy transports super- 
diffusively with a ~ 1.2. This is slower than ballistic trans- 
port ( a = 2 ) but faster than normal diffusion (a = 1 ). 
With temperature increases, the anharmonic terms appear due 
to excitation of the transverse vibrational mode. According to 
a theory proposed by Wang and Li [53], it is known that for a 
quasi- ID lattice model, the interaction between the transverse 
modes and the horizontal modes will result in a superdiffusive 
phonon transport. Moreover, Li and Wang [54] demonstrate 
a connection that the anomalous energy diffusion induces the 
divergent thermal conductivity. Combine the observed super- 
diffusion with a « 1.2 and their theory, the thermal conduc- 
tivity diverges with tube length for the tube length up to few 
micrometers, namely, the thermal conduction of SWNT does 
not obey Fourier's law. This provides the physical mechanism 
for the experimental and theoretical observed length depen- 
dent thermal conductivity of CNTs. 

In addition to CNTs, silicon nanowire (SiNW) is an- 
other promising one-dimensional nanomaterial that pushes the 
miniaturization of microelectronics towards a new level and 
have drawn significant attentions because of the ideal inter- 
face compatibility with the conventional Si-based technol- 
ogy. With the progressive research in its applications [55- 
57], more and more theoretical efforts have been made to 
understand the thermal properties of SiNWs due to the po- 
tential thermoelectric applications in both power generation 
and on-chip cooling [20, 58-60]. The impacts of temperature 
[61], surfaces structure [62-64], diameter [65, 66], tubular and 
core-shell structure [24, 67-70] and doping [71-73] have been 
reported. 

Like in CNTs, the length dependent thermal conductivity 
is also an important physical question to SiNWs. By using 
NEMD simulation, Yang et al. have studied the length de- 
pendence of thermal conductivity of SiNWs [9]. In order to 
establish a temperature gradient along the SiNW, the atoms 
close to the two ends are put into heat bathes with high and 
low temperature. Nose-Hoover and Langevin heat bathes [74] 
are applied to ensure the conclusions are independent of heat 
bath. The dependence of thermal conductivity on the length of 
SiNW is shown in Fig. 8. Both types of heat baths give rise to 
the same results. It is obvious that the thermal conductivity in- 
creases with the length as, k oc , even when the wire length 
is as long as 1.1/jm. Traditionally, the phonon mean free path 
(A) is a characteristic length scale beyond which phonons scat- 
ter and lose their phase coherence. In three-dimensional sys- 
tems, it gives Fourier's law when system scale L > A. By 
using the phonon relaxation time (~ lOps) in SiNWs, and 
the group velocity of phonon as 6400m/s, the mean free path 
A is about 60nm [9]. The maximum SiNW length (1.1/im) 
in this study is obviously much longer than the phonon mean 
free path of SiNW. This demonstrates that in SiNWs, Fourier's 
law is broken even when the length is obviously longer than 
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FIG. 8. (color online). The dependence of thermal conductivity of 
SiNWs on the longitude length. The black solid curves are the power 
law fitting curves (linear in log-log scale). For further details see in 
Ref. [9]. 



the traditionally mean free path. 

More interestingly, it is found that the length dependence of 
thermal conductivity is different in different length regimes. 
At room temperature, when SiNW length is less than about 
60nTO, the thermal conductivity increases with the length lin- 
early ( j3 « 1 ). For the NW with length larger than 60?ito, 
the divergent exponent j3 reduces to about 0.27. This criti- 
cal length (60nm) is just the value of mean free path, then 
this length dependent divergent exponent can be fully under- 
stood as below. There is weak interaction among phonons 
when the length of SiNW is shorter than mean free path, thus 
the phonons transport ballistically, like in the harmonic lattice. 
However, when the length of SiNW is longer than the mean 
free path, phonon-phonon scattering dominates the phonon 
transport and the phonon cannot flow ballistically. 

In addition, the diverged exponent j3 also depends on tem- 
perature. At 1000K, (3 is only about 0.15 when NW length is 
longer than 60nm. The decrease of j3 can be understood from 
the temperature dependent phonon interaction. In SiNWs at 
high temperature, the displacement of atoms increases, which 
induces more phonon-phonon interaction, and reduction in j3. 

The divergence of thermal conductivity can be understood 
from the phonon density of states spectra (PDOS) as shown 
in Fig. 9. In the left column of Fig. 9, the higher energy part 
(short wavelength) in PDOS spectra is not sensitive to the wire 
length. The peak at about 16 THz in PDOS of bulk silicon ap- 
peares even when the NW length is only 4 nm. However, as 
shown in right column of Fig. 9, the lower frequency part in 
PDOS spectra is much different from that of bulk Si. In con- 
trast to the continuous spectra of bulk silicon, there are many 
discrete peaks in PDOS spectra of NW. For a short SiNW, 
the low energy phonon density is very low, thus the thermal 
conductivity is low. With the increase of length, more and 
more (long wavelength) phonons are excited, which leads to 
the increase of thermal conductivity. However, due to the size 
confinement effect, the energy density of acoustic phonon in 
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FIG. 9. (color online). The phonon density of states along the lon- 
gitude direction of SiNWs with different lengths, and the phonon 
density of states of bulk Si. For further details see in Ref. [9]. 
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FIG. 10. (color online). Energy diffusion in SiNW at room tempera- 
ture and at 1000K. The length of SiNW is 140 nm. For further details 
see in Ref. [9]. 



PDOS spectra of SiNW is much smaller than that of bulk Si, 
and leads to the low thermal conductivity of SiNW compared 
to that of bulk Si. 

So far we show the length dependent thermal conductivity 
of SiNWs. Now let us turn to the energy diffusion process in 
SiNWs. The heat energy diffusion in SiNW was firstly sim- 
ulated by Yang et al. in 2010 [9]. Followed the approach 
applied in diffusion in CNT by Zhang and Li [8], they first 
thermalize the system to an equilibrium state with temperature 
To, then atoms in the middle layer are given a much higher 
temperature T\ . The evolution of the energy profile along the 
chain is then recorded afterwards. Fig. 10 shows (V 2 (i)) 
versus time in double logarithmic scale, so that the slope of 
the curve gives the value of a. For the NWs of length of 
140nm, we obtain a = 1.15 and 1.07 at 300K and 1000K, re- 
spectively. Different theory describes the physical connection 
between energy diffusion and thermal conductivity [54, 75]. 
For instance, in normal diffusion, the phonon transports dif- 
fusively, it corresponds to a size-independent thermal conduc- 
tivity. This is what we have in bulk material and the Fourier's 
law is valid. However, in the ballistic transport, the thermal 
conductivity of the system is infinite when the system goes 
to thermodynamic limit. This is the case for one-dimensional 
harmonic lattice. In another case, if a < 1, which we call 
sub-diffusion case, corresponds to j3 < 0, namely, the ther- 
mal conductivity of the system goes to zero. Thus the system 
is an insulator. In the super diffusion regime, it predicts that 
the thermal conductivity increases as the length of the system 
increases. This is what we observed in SiNWs. 

Thus the anomalous heat diffusion is responsible for the 
length dependent thermal conductivity in SiNWs. Combine 
with the results in CNTs, all these reports provide strong evi- 
dence that Fourier's law of heat conduction is not valid in one 
dimensional and quasi one dimentional nanostructures. More- 
over, these theoretical studies demonstrate that nanowire and 
nanotube are promising platforms to verify phonon transport 
mechanisms. 



Computational Study on Thermal Conductivity of Graphenes 

Besides to the one-dimensional nano materials, Graphenes 
[28, 76] have attracted immense interests recently, mostly 
because of their unusual one-atom-thick structure. Superior 
thermal conductivity (as high as 5000 W m~l K~l) has 
been observed in graphene [30], which has raised the excit- 
ing prospect of using them for thermal devices. 

Based on graphene, narrow (sub-10 nm) graphene nanorib- 
bon (GNR) can be obtained either by cutting mechani- 
cally exfoliated graphenes, by patterning epitaxially growth 
graphenes or by unzipping carbon nanotubes (CNTs) [28, 77, 
78]. GNRs are materials with distinctive electronic and ther- 
mal transport properties and are candidates for field effect 
transistors of future carbon-based nanoelectronics. There are 
rich physical phenomena about thermal property of GNRs. 
For instance, intrinsic anisotropy originate from different 
boundary condition at ribbon edges is observed in narrow 
GNRs, results to the room temperature thermal conductance 
of zigzag GNRs is about thirty percent larger than that of arm- 
chair GNRs [79, 80]. The effects of edge, roughness, and hy- 
drogen termination have been investigated by different simu- 
lation methods [81-86]. And isotopic doping can reduce re- 
markably the thermal conductivity of graphene [87]. In addi- 
tion to the isotope randomly doping, the impact of disorder de- 
fect on thermal conductance [88, 89] and thermoelectric figure 
of merit of GNRs [90] has been investigated by using density 
functional theory calculations combined with the nonequilib- 
rium Green's function method. It is found that the figure of 
merit can be remarkable enhanced five times by randomly in- 
troducing vacancies to the graphene. And remarkable ther- 
mal conductivity enhancement has been found in graphene 
nanoribbons under homogeneous uniaxial stretching due to 
a lot of dispersive phonon modes are converged to the low 
frequency region [91]. And thermal conductivity of GNRs 
is very sensitive to tensile strain [92]. More interesting, the 
presence of a substrate [93, 94] and layer-layer interaction 
[95-98] lower the thermal conductivity. The excellent thermal 
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property of GNRs makes them ideal candidate for nanoscale 
phononic devices [95, 99-101]. 

Similar to the cases in CNT and SiNW, thermal conduc- 
tivity of GNRs also depends on their sizes. By using equi- 
librium molecular dynamics simulations, the effect of GNR 
width on thermal conductivity has been studied [81]. They 
found that the thermal conductivity increases with width, but 
for width larger than hnm the thermal conductivity becomes 
relatively size independent. The larger conductivity values 
at very small sizes is explained from the limited number of 
phonons in the system. In their simulation, the width range is 
from 2nm to Ylnm. The similar width dependence of thermal 
conductivity was also observed by non-equilibrium molecular 
dynamics simulations [102]. Moreover, thermal conductivity 
also depends on GNR length. In their work, for GNR with 
width of 2nm, its thermal conductivity increases from about 
200 W/mK to 400W/mK when length increases from lOnm to 
60nm. Based on the extrapolation procedure, they predicted 
the thermal conductivity of GNRs with a length of 2 microme- 
ter will be about 3000W/mK, which is close to the experimen- 
tal reports. However, as the maximum length in their simula- 
tion [102] is only about 60?ito, it is not long enough to pro- 
vide detailed information about the divergent exponent. Fur- 
ther simulation studies on much larger GNRs will be greatly 
helpful to explore whether Fourier's law is still validity in two- 
dimensional one-atom-thick systems (See also [103]). 



HEAT TRANSPORT IN LOW DIMENSIONAL MODELS 



The anomalous heat transport phenomena observed in the 
above mentioned low dimensional nanostructures can be par- 
tially understood by general theories on low dimensional 
models, which generally demonstrate that the phenomeno- 
logical Fourier's law may not be valid. Actually, ever since 
Fourier proposed this law on heat conduction about 200 year 
ago, enormous number of studies has been conducted to 
fully understand it from fundamental physics. But a rigor- 
ous derivation from microscopic Hamiltonian dynamics is still 
absent. However, from these studies, we are already able to 
draw the surprising conclusion that heat conduction in low di- 
mensional systems may not follow the Fourier's law. In this 
section, we will briefly review the fundamental studies on heat 
transport in low dimensional systems: from theories to numer- 
ical simulations on the toy models. Especially, we introduce 
the works that study anomalous heat conduction from the as- 
pect of anomalous energy diffusion. 

In order to understand the realistic models from a rigor- 
ous perspective, it is better to begin with the simplest mod- 
els which involve only the basic but important ingredients, 
such as linearity and disorder, that will possibly lead to the 
expected normal or anomalous transport behaviors. 



Low Dimensional Harmonic Lattices 

The ordered harmonic lattice is the simplest lattice model 
that we can start with. From the early Boltzmann transport 
equation, it is not surprising that this linear model does not 
have well defined transport coefficients due to the lack of in- 
teraction between phonon modes. Therefore, heat transport in 
harmonic models should not obey Fourier's law. Actually, a 
clearer clue was found after Debye extended Boltzmann ki- 
netic theory of ideal gas and established the expression of 
heat conductivity k = cvl, where c is the heat capacity, v the 
phonon velocity, and I the phonon mean free path, because the 
non-interacting phonons should have infinite mean free paths. 

The first explicit results on heat transport in classic har- 
monic models were given by Rieder, Lebowitz and Lieb 
[1]. They studied a harmonic chain connected to stochas- 
tic Langevin heat baths. The lattice has Hamiltonian H = 

Y^iLi 2m + Szl^ 1 2~k{xi+\ — xi) 2 . The particles at the ends 
are subjected to additional forces F^jv = r)i,N — 7Ui,jv, 
where the noise terms rji t N are independent Gaussian random 
processes with mean zero and variances (171 jv(t)f/i,jv(0)) = 
2-fkBTi^S(t). Rieder et al. proved that, at large N, the heat 
flux will saturate to [1] 



"27" 



(n-T N ), 



(6) 



where v = This heat flux is independent of the length TV 
for large N, which demonstrates the expected ballistic trans- 
port. Another important result on harmonic model is that the 
non-equilibrium temperature profile is flat at the center part 
with temperature T = (T\ + Tjv)/2. Temperature jumps 
occur at the boundaries. In this sense, temperature gradi- 
ent cannot be established in harmonic lattices. Another work 
by Nakazawa extended Rieder et aZ.'s model by introducing 
on-site harmonic potential to all sites and a similar length- 
independent heat flux was found [2]. In the same work [2], 
Nakazawa also proved that high dimensional harmonic lat- 
tices can be reduced to a ID problem. Therefore the con- 
clusions are similar. 

Due to the absence of scattering mechanism in harmonic 
models, heat transport in these lattices is ballistic. Therefore, 
it is necessary to introduce scattering of phonons to get the de- 
sired diffusive transport. Basically, two approaches are com- 
monly adopted: by introducing disorder or nonlinearity to the 
system. Both of these two methods can cause phonon-phonon 
interaction. However, as we will see later, neither method 
works effectively for diffusive transport in low dimensional 
lattices. 

For the first approach, we consider the harmonic lattices 
with disorder. Disorder can be generally introduced by ran- 
domly assigning varies particle masses or spring constants, or 
both. These two types of disorder do not change the underly- 
ing physics very much. The former with random masses is just 
like isotopic doping, which is more interesting and relevant in 
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applications. Therefore, it is studied intensively. The pres- 
ence of disorder will generally cause localization of normal 
modes, which could be partially understood as an analogue to 
the Anderson localization of electrons. In the Anderson tight- 
binding model, all the eigenstates of electrons are localized 
in one and two dimensional cases, therefore the systems are 
electric insulators. However, in the phonon case, the picture 
is much more complicated. 

For one dimensional disordered harmonic chains, Allen and 
Ford first observed that the thermal conductivity of an infi- 
nite disordered harmonic chain is finite [ 104], unlike the An- 
derson model of electrons in which the conductivity will ex- 
ponentially decay to zero. Furthermore, the eigenstates of 
electrons in the strongly disordered cases are all exponen- 
tially localized. However, in disordered media, Matsuda and 
Ishii [105] showed that only the high frequency phonons has 
this property, while phonons with low frequencies oj < u>d 
can be considered as delocalized phonons, since the localiza- 
tion lengths are greater than the system length. The charac- 
teristic frequency Ud depends on the variance of mass dis- 
tribution (Am 2 ) as Lu d ~ ( jv(A™t 2 ) ) 1 ^ 2 ' wnere ^> m an d 
N are the spring constant, the average mass and the lattice 
length, respectively. Therefore, low frequency phonons will 
contribute to heat transport, which prevents the disordered 
harmonic lattices from being a thermal insulator. Matsuda 
and Ishii also studied the thermal conductivity of disordered 
chains connected to white noise heat baths (Langevin baths) 
and baths modeled by ordered semi-infinite harmonic chains 
(Rubin baths), respectively [105]. They found that for both 
models (a): disordered harmonic chain connected to Langevin 
baths with fixed boundary conditions; and model (b): dis- 
ordered harmonic chain connected to Rubin baths with free 
boundary conditions, their thermal conductivity will diverge 
with the system length as n ~ with /? = 1/2. However, 
Casher and Lebowitz later rigorously proved that the correct 
exponent for model (a) should be /3 = —1/2 [106], while for 
model (b) the result /3 = 1/2 is supported by numerical simu- 
lations [107] and a rigorous proof was provided by Verheggen 
later [108]. 

Altogether, the thermal transport properties of disordered 
harmonic chains would probably depend on the both the 
boundary conditions and the heat baths. To clarify this point, 
Dhar restudied the thermal transport in disordered harmonic 
chains connected to various baths modeled by generalized 
Langevin equations [109]. Using a Langevin equations and 
Green's function (LEGF) formalism, Dhar found that in disor- 
dered harmonic chains, the exponent j3 depends on the spec- 
tral properties of the baths. The previous model (a) and (b) 
are only two special cases of the generalized Langevin heat 
baths, which have j3 = —1/2 and 1/2, respectively. A spe- 
cial choice of heat bath spectrum can even lead to normal 
heat transport with (3 = 1, obeying Fourier's law [109]. A 
more detailed study on the effects of the boundary conditions 
on the transport properties were carried out by Roy and Dhar 
later using the same LEGF formalism [110]. The results are 
a little contradictory to those of [109], which states that the 




100 1000 10000 100000 

N 



FIG. 11. (J) vs. N for different boundary conditions with tempera- 
ture bias fixed. Results are given for the Langevin heat baths (model 
a) and Rubin baths (model b). The two straight lines correspond to 
the asymptotic power law behaviors J ~ TV -1 ' 2 and J ~ iV~ 3/ ' 2 , 
respectively. For further details see Ref. [110]. 

different exponents /3 are indeed dependent on the different 
boundary conditions, instead of the spectral properties of the 
baths. Clear numerical evidence also supports this conclusion 
as shown in Fig. (11), i.e., for disordered harmonic chains 
with free boundary conditions j3 = 1/2 and with fixed bound- 
ary conditions j3 = —1/2, no matter whether Langevin baths 
or Rubin baths are used. Ref. [110] also showed that, in the 
low frequency regime, the transmission coefficient of the dis- 
ordered chains T(oS) can be approximated by the transmission 
of the ordered chains /"ordered (k>). Therefore, heat flux in a 
disordered harmonic chain can be estimated as 

J ~ (T L - T R ) / Ordered (<*>) dui, (7) 

Jo 

where LUd is the low frequency boundary ujd ~ ( N(/\m*) Y^ 2 - 
Considering the fact that the transmission coefficients for low 
frequency modes scale as T or d e red(^) ~ 1 for free boundary 
conditions and T or dered(w) ~ ^> 2 for fixed boundary condi- 
tions, respectively, we can obtain Jf ree oc (Tj, — Tr) ojd oc 
(Tj, — Tjij/N 1 / 2 for free boundary conditions and Jfi XG d oc 
{Tl - T R ) cuj oc (T L - T R )/N 3 / 2 for fixed boundary condi- 
tions, respectively. 

Till now, we only discussed the disordered harmonic chain 
without onsite pinning potential (except the case of fixed 
boundary conditions, the particles at the ends are pinned). The 
absence of pinning in this model makes the system transla- 
tional invariant in a coarse grained manner. As a consequence, 
the phonon modes are not fully localized and the low fre- 
quency ones can transport heat. When onsite potentials are 
involved, the picture is a little different. In their work [111], 
Dhar and Lebowitz proved that with onsite potentials, all the 
phonon modes are localized. Consequently, the heat flux J 
decays exponentially with the chain length and the system is a 
thermal insulator in the thermal dynamic limit. This could be 
understood because the onsite potentials will typically break 
the translational invariance of the chain and eliminate the low 
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frequency phonon modes in this structure. The remaining 
phonon modes have high frequencies and are localized. 

Compared with the one dimensional lattices, the studies on 
two dimensional disordered ones are much less due to the 
mathematical difficulties. In a renormalization group study by 
John et al. [112], it was found that the low frequency modes 
in 2D disordered harmonic lattices are not localized either 
as in ID cases. The localization length diverges as e 1 '^ in 
the low frequency regime. Therefore, modes with frequency 
uj < uid = \og ' 2 {N) can be considered as extended, which 
will contribute to heat transport. 

In another study by Lee and Dhar[113],2D disordered har- 
monic lattices connected to two types of stochastic heat baths 
are investigated. The first type is modeled by uncorrelated 
Gaussian processes and the other by exponentially correlated 
Gaussian processes. In their work, disorder is introduced as 
isotopic doping. A mass M = 2 is randomly assigned to half 
of the particles, while the rest have mass m = 1. The authors 
simulated square lattices with widths up to L = 256 and found 
a power law divergent heat conductivity in both two types of 
models. For uncorrelated baths the exponent /3 « 0.41 while 
for correlated baths (3 « 0.49. In the same work [113], Lee 
and Dhar also studied a special case of correlated disorder in 
which the lattice is only disordered in the conducting direc- 
tion. It was found that this model can be transformed into 
an effective one dimensional problem thus being mathemat- 
ically tractable. Analytical results showed that the exponent 
j3 = —1/2 which was also verified by numerical simulations 
using square samples with widths up to L = 512. In another 
work [114], Yang considered 2D lattices with bond-missing 
defects connected to Nose-Hoover heat baths. The results 
showed that when the defect density is large enough, temper- 
ature gradient can build up and the heat conductivity is finite, 
while for the case of small defect density, the conductivity will 
diverge logarithmically. 



Low Dimensional Anharmonic Systems 

It is clear now that disorder alone cannot effectively make 
the phonon transport in the harmonic chains diffusively. As 
an alternative approach to introduce phonon-phonon scatter- 
ing, the effect of nonlinear interaction on the properties of 
heat transport in low dimensional systems has also been inten- 
sively studied, especially by numerical simulations. Because 
of the intrinsic non-integrability of the system Hamiltonian, 
analytical results are rare. 

Basically, there are three different approaches to study heat 
transport analytically: the mode coupling theory (MCT) [4, 
115], the renormalization group theory [45] and the Peierls- 
Boltzmann kinetic theory [116]. All are based on Green-Kubo 
formula [117, 118] 



total heat flux J. 

Strictly speaking, Green-Kubo formula can only be applied 
to infinite systems, as indicated by the limit taken lirnx_>oo, 
which should be performed before lim^oo strictly. However, 
as we already know, in many low dimensional systems, the 
heat conductivity in thermal dynamic limit is divergent. The 
focus of such anomalous situations is how the conductivity 
diverges with the system size. To utilize the Green-Kubo for- 
mula to study the length dependence of heat conductivity in 
these cases, a usual procedure is to truncate the upper limit of 
the integration in Eq. (8) to t c ~ L/v s . This is commonly be- 
lieved as that the sound waves propagate to the boundaries at 
a finite speed v s , which will lead to a fast decay of correlation 
C(t) at time ~ L/v s [14, 15]. Using this treatment, the pre- 
viously mentioned three theories, i.e. mode coupling theory, 
renormalization group theory and Peierls-Boltzmann equation 
kinetic theory, are all aimed to calculate the long time tail of 
the correlation function C(t). A basic result shared by all 
these theories is that the systems with nonlinear interaction 
remain anomalous except those subjected to onsite potentials. 
For ID cases without onsite potential, these theories predict 
that the conductivity of the system would diverge as a power 
law k ~ N 13 . However, the values of (3 differ from one to an- 
other. For 2D cases, the mode coupling theory and the renor- 
malization group theory predict logarithmic divergence. 

The mode coupling theory approach was first used by Lepri 
et al. [4, 115]. The basic idea is that the divergence of heat 
conductivity is due to the slow relaxation of phonon modes 
with long wavelengths. A later development of this theory 
was done by Delfini et al. [119]. By constructing a correlation 
function G(k,t) = (Q*(k, t)Q(k, 0)) of the normal mode 
coordinates Q(k) = J2 n =i Xn exp(— ifcn), the authors 
were able to get the exact evolution equations for G(k, t) us- 
ing the Mori-Zwanzig projection approach [115]. Using some 
approximations, these equations were solved self-consistently 
for the Fermi-Pasta-Ulam (FPU) chains with interaction po- 
tential U(x) = k 2 x 2 /2 + k 3 x 3 /3 + fc4£ 4 /4, which gives 
G(k,t) = A(k,t)e lul< - k ' ,t + c.c. for a small wavenumber k. 
In this formula, uj(k) is the temperature dependent disper- 
sion relation which is obtained using harmonic approximation 
Lu(k) = 2| sin(fc/2)|, and A(k, t) has the form 



A{k,t) 



g{y/3k 2 k B T/2Trtk 3 / 2 



k 3 + 0; 



— L- lim lim / C(t)dr, 

ksl t-¥0O L^roo Jq 



(8) 



where C(t) = (J(t)J(0)) /L is the time correlation of the 



g{ v / 15(k 4 k B T/2n)' 2 tk 2 ) k 3 = and k 4 ^ 0. 

(9) 

Finally, the heat flux correlation function C (t) was obtained 
fromC(i) oc Efc(^sr) 2G2 (M)- It was shown that C(t) ~ 
t" 2 / 3 for the case k 3 ^ and C(t) ~ t' 1 / 2 for k 3 = but 
&4 7^ 0. Therefore, by inserting C(t) into the Green-Kubo for- 
mula and adopting the cut-off time t oc N, we obtain (3=1/3 
and (3 = 1/2 for each case, respectively. The mode coupling 
theory can also provide predictions for higher dimensional 
systems. It was shown that [14] for two dimensional systems, 
C(t) ~ t -1 , so that the heat conductivity would diverge log- 
arithmically with system size. While for three dimensional 
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systems, C(t) ~ t~ 3 / 2 , with which we obtain normal heat 
conduction. 

In a later work, Wang and Li investigated the effect of trans- 
verse degrees of freedom in one dimensional chains using this 
mode coupling theory [53]. The system under investigation 
has Hamiltonian 

H = E + ^(i r «+i - r 'i - °) 2 + c 10 ) 

where (r^, r^) is the canonical coordinates, and <fri is the angle 
between vectors r^+i — and r,_i — r^. This model can be 
regarded as a simplification of more realistic polymer chains. 
The model is still one dimensional but the motions of particles 
are two dimensional. Both numerical simulations and mode 
coupling analyses suggest that in the presence of interaction 
with transverse modes, j3 = 1/3. When the coupling between 
transverse modes and longitudinal modes is weak, j3 = 2/5. 

The second approach was proposed by Narayan and Ra- 
maswamy which uses hydrodynamic equations and renormal- 
ization group theory [45]. It was argued that in an interacting 
ID system with large system length, the thermal fluctuation 
will wipe out the long-range order, and consequently the sys- 
tem will behave like fluid. Therefore it is expected that this 
approach can be used generally to describe heat transport in 
ID phenomenon. By assuming that the only conserved quan- 
tities in the system are the total number of particles, the to- 
tal momentum and the total energy, one can obtain three hy- 
drodynamic equations describing the evolution of the particle 
density field and the velocity field, with addition terms de- 
scribing the thermal noise. Narayan and Ramaswamy then 
solved these equations using linear response approximation 
and finally obtained the thermal flux correlation by consider- 
ing the symmetries of the system. It was demonstrated that 
the flux correlation function C(t) ~ i~ 2 / 3 asymptotically for 
the ID case. Therefore, one obtains /? = 1/3 by the cut-off 
time reasoning. For 2D systems, the authors showed that the 
conductivity would diverge logarithmically with system sizes 
[45]. 

In another work by Mai and Narayan [120], the authors 
claimed that the renormalization group theory analyses can 
be applied to lattice models as well, even for stiff one dimen- 
sional oscillator chains including the FPU model. The con- 
ductivity of ID lattices should also diverges in the same way 
as fluid, i.e. diverges in a power-law with exponent (3 = 1/3. 
However, in a consequent work by Hurtado [121], the break- 
down of hydrodynamics in a simple one dimensional fluid was 
reported. 

The estimation of the exponent /3 in anomalous heat trans- 
port using Peierls-Boltzmann equation is first carried out by 
Pereverzev [116]. The Boltzmann equation is originally used 
to describe the the phase space density evolution of kinetic 
gases. From it we have already known that the heat conduc- 
tivity can be written as k = cvl, where c, v and I are the 
specific heat, sound velocity and mean free path of the gases, 
respectively. Peierls developed this theory to describe the 



phonon transport in solids and found a similar formula for the 
conductivity k ~ J CfcU 2 Tfcdfc, where k specifies the phonon 
modes and Tfc is the phonon relaxation time. Using this 
Peierls-Boltzmann approach, Pereverzev studied the FPU-/? 
mode [116]. Considering that anomalous heat conduction is 
mainly due to the ballistic-like behavior of phonons with small 
wavenumber k. the dependence of on the wavenumber k 
in this regime is the main interest. After some approximation, 
Pereverzev was able to obtain ~ fc~ 5 / 3 for small k. It was 
also found that heat flux correlation function can be calculated 
by C(t) = 2k] ^ T Jq e~ t / Tk vfdk. Because Vh is almost con- 
stant for small k, one obtains C(t) ~ i~ 3 / 5 and /3 = 2/5. 
This result was confirmed later by Lukkarinen and Spohn us- 
ing a more rigorous Peierls-Boltzmann approach [122]. For 
momentum non-conserving cases, a finite heat conductivity 
was obtained by Aoki et al. [123] using this approach. 

Compared to the few analytical results on low dimensional 
heat conduction, numerical results are much more abundance. 
Here we only briefly highlight the results on some well studied 
lattice oscillator models. 

Due to the landmark work by Fermi, Pasta and Ulam, 
the FPU model has received great attraction in understand- 
ing nonlinear statistical mechanics. In the FPU-/? model, 
interaction between nearest particles has the form U(x) = 
k2X 2 /2 + k4X 4 /A. This model is one of the simplest lattice 
models which can be used to investigate the effect of nonlin- 
earity on heat transport. 

The first numerical study on heat conduction of FPU-/? 
model was carried out by Kaburaki and Machida [124]. In 
their study, the ends of FPU-/? lattices were connected to 
Boltzmann heat baths at different temperatures. Detailed cal- 
culations showed that linear temperature profiles can be es- 
tablished in the steady state. Moreover, it was found that the 
conductivity tends to saturate when the chain length increases. 
However, it is clear now that neither of these findings reveals 
the true heat transport behaviors of FPU-/? lattices. 

The first indication of anomalous heat conduction violat- 
ing Fourier's law in FPU-/? lattices was presented by Lepri 
et al. [3]. FPU-/? lattices connected to Nose-Hoover heat 
baths with lengths up to N = 400 was studied. It was found 
that the heat conductivity diverges as k ~ N 13 with exponent 
/? = 0.55 ± 0.55. Another finding in this work is that the tem- 
perature profile is non-linear even when the temperature bias 
is very small, a result contradictory to the prediction from the 
Fourier's law. In a subsequent work [4], the divergent expo- 
nent was improved to /? = 0.37 with length up to N = 2048. 
These pioneering works in one dimensional lattices have in- 
spired enormous subsequent studies on anomalous heat con- 
duction in low dimensional systems. 

It should be mentioned that the system lengths in the studies 
above are too small to draw any quantitative conclusions. Ac- 
tually, in non-equilibrium simulations, the heat conductance 
will be greatly affected by the contact resistances between the 
lattice and heat baths. Therefore, it is necessary to go to a 
very large system length to make the contact resistances neg- 
ligible. In a recent work [125], Mai et al. have carefully taken 
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FIG. 12. Heat conductivity kne vs. lattice length N in log-log scale. 
In the rightmost part of the figure, it can be seen that the running 
slope increases as the system length is increased. Solid and dashed 
lines with slope 2/5 and 1/3, respectively, are drawn for reference. 
For more details please see Ref. [126]. 



this into consideration. They studied the FPU-/? chain with 
length up to N = 65536 using two different types of heat 
bath: the stochastic Langevin heat baths and the determinis- 
tic Nose-Hoover heat baths. It was found that when the sys- 
tem size is small, the calculated thermal conductances differ 
a lot when different heat baths are used. Meanwhile, when 
the system size is large enough, the discrepancy between the 
calculated heat conductances is negligible, which ensures us 
that the boundary effects could be omitted. In this situation, it 
was then found that the exponent j3 decreases to a final result 
j3 = 0.333 ± 0.004. The authors claimed that this result is 
consistent with the renormalization group theory [45] which 
supports its universality. 

However, in a quite recent work, Wang et al. [126] restud- 
ied the FPU-/3 model and the purely quartic model with an 
even longer length up to N = 131072. The interaction in 
purely quartic model is U(x) = \k^. This model is the 
high temperature limit of the FPU-/3 model. They found that 
the slope dln/c/dlniV is not monotonically decreasing as 
shown by Mai et al. [125]. The slope starts to increase when 
N is close to 10 5 and finally reaches j3 = 2/5 instead of 1/3 
(see Fig. 12). In order to remove the boundary effects, the 
authors also performed equilibrium simulations without heat 
baths to calculate the heat flux correlation function C(t) in 
the Green-Kubo formula. Again, the slope d In C(t)/d In N 
increases when t is very large. The heat flux correlation fi- 
nally saturates at C(t) ~ i~ 3 / 5 , which supports their non- 
equilibrium results (3 — 2/5 and the Peierls-Boltzmann theory 
[116]. 

For one-dimensional lattices with momentum conserving 
interparticle interactions, all theories predict a divergent heat 
conductivity n oc L 13 with universal j3. However, in a most 
recent study on FPU- (3 lattices with nearest-neighbor (NN) 
and next-nearest-neighbor (NNN) coupling [127], Xiong et al. 
found that the exponent (3 strongly depends on the ratio of the 
NNN coupling to the NN coupling 7. The Hamiltonian of the 
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FIG. 13. Dependence of ,9 on 7 in the model with Hamiltonian Eq. 
(11). For more details please see Ref. [127], 



system under investigation is 
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where V(x) is of FPU-/? type ^x 2 + \x . In their study, the 
authors used the reverse nonequilibrium molecular dynamics 
method (RNEMD) [128] to simulate the system. Firstly, this 
method was tested under condition 7 = and length L = 
2496 which showed that the temperature gradient was well 
established and f3 = 0.325 ± 0.002. Consequent simulations 
using different 7 revealed that the exponent f3 varies continu- 
ously with 7. Starting from f3 = 0.325 when 7 = 0, (3 keeps 
decreasing until reaching its minimum 0.25 at 7 = 0.25. Af- 
ter this point, (3 increases to 0.35 at 7 = 1 as shown in Fig. 13. 
Therefore, it was suggested that a universal exponent (3 does 
not exist. However, it should be mentioned the result /3 at 
7 = is quite different from Wang et al. 's which was calcu- 
lated using both non-equilibrium method with Langevin heat 
baths and equilibrium method with Green-Kubo formula with 
length up to 131072 [126]. 

For momentum non-conserving models, it has been men- 
tioned that all theories predict a finite heat conductivity. Nu- 
merical simulations also support this conclusion. Actually, 
before the three theoretical predictions were announced, Hu 
et al. had realized that momentum non-conservation should 
be a necessary ingredient to obtain the Fourier's law from their 
work on the Frenkel-Kontorva (FK) model with a periodic on- 
site potential V(x) oc cos(aa;) in addition to the harmonic in- 
terpartical interaction [129], Two other studies, one by Aoki 
et al. [130], and the other by Hu et al. [131], on </> 4 model 
with onsite potential V(x) = k^x^/A further confirms this 
conclusion. 

What is the disorder effect to anharmonic lattices? This 
problem was first investigated by Pay ton et al. [132]. They 
performed non-equilibrium simulations on disordered anhar- 
monic lattice connected to stochastic baths. It was found that 
the anharmonicity could greatly enhance the heat current. But 
due to the limited computer facility at that time, they were not 
able to verify the Fourier's law. The first systematic study of 
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the disorder effect in anharmonic lattice was taken by Li et al. 
[133]. In their work, Li et al. studied the FPU chain with 
mass disorder connected to Nose-Hoover heat baths. Their 
results showed that at low temperatures, the system obeys the 
Fourier's law. While at high temperatures, the conductivity 
diverges with exponent j3 = 0.43. However, in a later study 
by Dhar and Saito [134], the authors found that the exponent 
j3 is still 1/3 at low temperatures if one goes to much larger 
system size. Therefore, there is no crossover from diffusive to 
superdiffusive in disordered FPU chain. Disorder plays less 
important role in j3 in anharmonic chains. 

For disordered anharmonic chains with pinning, Dhar and 
Lebowitz studied the general Hamiltonian systems [111] 
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in which xq = xjsr+i = 0. It was found that when anhar- 
monicity is absent, which corresponds to the disordered har- 
monic chain with harmonic onsite potential, the heat flux de- 
cays exponentially with length and the system is a thermal 
insulator. However, introduction of a small will lead to the 
J ~ l/N dependence, resulting in a diffusive heat transport 
arises. 

Regarding the two dimensional oscillator lattices, only few 
numerical studies are carried out to investigate the size depen- 
dent heat conductivity. But the good news is that predictions 
from theories do not contradict one another. Both the renor- 
malization group theory and the mode coupling theory predict 
a logarithmic divergence for the heat conductivity. 

For numerical investigations, Lippi and Livi studied two di- 
mensional lattices with size N x x N y [135]. Heat was trans- 
ported along x direction. The Hamiltonian of the system has 
the form 
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where ,xy are the canonical coordinates 



(13) 

The interaction 

was taken to be the FPU-/3 type U (x) = x 2 /2 + k^x 4 /^ and 
the Lennard-Jones type U(x) = A/x 12 - B/x e + B 2 /4A, 
respectively. The parameters in the Lennard-Jones potential 
were chosen so that its Taylor expansion at the minimum coin- 
cide with the FPU-/3 potential. Using non-equilibrium method 
with Nose-Hoover heat baths, Lippi and Livi firstly demon- 
strated that with increasing Ny and fixed N x , the current sat- 
urates at a small ratio N y /N x . Subsequent calculations were 
performed using fixed ratio N y /N x = 1/2. It was found that 
the conductivity will diverge logarithmically with N x for both 
models. Additional equilibrium simulations showed that the 
heat flux correlation functions decay as t~ x . Using Green- 
Kubo formula, we again obtain n ~ ln(iV x ). This result is 
consistent with predictions from theories. 



There also exists numeric works that contradicts the predic- 
tion of logarithmic divergence. Yang and Grassberger stud- 
ied a 2D Hamiltonian system involving FPU type interaction 
[136]. The Hamiltonian of the systemis similar to Eq. (13) but 
with scalar displacements. The results showed that when the 
ratio N x /N y is small which corresponds to the 2D behavior, 
the conductivity diverges exponentially with (3 = 0.22 ±0.03. 
When N y = 1 which is actually a ID system, /3 = 0.37±0.01. 
In another recent work, Shiba and Ito also studied the same 
system as Eq. (13) with the FPU-type interaction [137]. How- 
ever, a power law divergence was obtained by using lattices 
with sizes up to 384 x 768. It is found the conductivity also 
diverges in a power-law with /3 « 0.268, even it was assured 
that the boundary effect is negligible in that size. 



Anomalous Energy Diffusion in Low Dimensions 

The nonlinearity of the interaction in the lattices prohibits 
us to analytically solve most of the models. Instead, a lot of 
quasi-one-dimensional gas channel systems consists of non- 
interacting particles have been proposed to study heat trans- 
port in low dimensions. These works helped a lot in under- 
standing the basic ingredients required by the Fourier's law. 
Moreover, they helped to understand heat transport from an- 
other aspects of energy diffusion. 

Alonso et al. first studied the heat conduction in Lorentz 
gas channels [138]. The model consists of a quasi-lD billiard 
with periodically distributed semicircular scatterers as shown 
in Fig. 14(a). Particles carrying heat move along the billiard 
and exchange energy with Maxwell thermal reservoirs at the 
ends. In this model, no particle can move between the two 
reservoirs without being scattered by the semicircles. Due to 
the exponential separation of trajectories, this system has a 
positive Lyanpunov exponent, i.e., it is chaotic. Alonso et al. 
verified that the heat conduction obeys Fourier's law which 
was explained by deterministic particle diffusion with energy 
dependent diffusivity D[E) — E 1/2 . 

Later, in order to resolve the role of dynamic chaos on de- 
termining the heat transport behavior, Li et al. studied the 
Ehrenfest gas channels with isosceles right angles replacing 
the semicircles in the Lorentz gas channels [139]. Their con- 
figurations are shown in Figs. 14(a) and 14(b). The linear 
separation of nearby trajectories in these models make them 
different from the Lorentz gas channels in underlying dynam- 
ics. It was found that when the scatterers are periodic, the heat 
conductivity diverges with system length as k ~ N om . Nor- 
mal heat conduction k ~ N° can only be reached when either 
position disorder (Fig. 14(b)) or size disorder (Fig. 14(c)) is 
introduced. These results imply that chaos is not necessary for 
normal heat conduction and Fourier's law. Moreover, the dif- 
fusive behavior of the particles are also investigated using the 
mean square displacement. Numerical simulations showed 
that for the disordered case, the mean square displacement in- 
crease linearly with time, while the periodic case is superdif- 
fusive (Ax{t) 2 ) ~ t a with a = 1.672. 
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FIG. 14. The geometrical configuration of billiard gas channels, 
(a) The Lorentz gas channel, (b) The position-disordered Ehrenfest 
channel, (c) The height-disordered Ehrenfest channel, (d) The polyg- 
onal gas channel, (e) The triangle gas channel. 



Subsequently, two additional works on quasi- ID billiard 
model but with different scatterers are studied by Alonso 
et al. [140] and by Li et al. [141], respectively. In the 
first one, the scatterers have polygonal shape (Fig. 14d). By 
fixing <pi = (V5 — l)7r/8 and adjusting <p2 = n/q with 
q = 3,4, 5, 6, 7, 8, 9, Alonso et al. found that the time de- 
pendence of the MSD (Ax(i) 2 ) ~ t 13 is superdiffusive for 
q = 3 ((3 = 1.30), subdiffusive for q = 4 (,5=0.86) and nor- 
mal for 5 = 5,6, 7, 8, 9 (/3 ~ 1). Numerical simulations also 
show that for q = 5,6,7,8,9, the heat conductivity of the 
system converges to a non-zero constant in thermal dynami- 
cal limit, while it diverges for q = 3 and decays to zero for 
q = 4 [140]. In the second work [141], the scatterers are 
triangles (Fig. 14). It was found that the transport behaviors 
of this model depend sensitively on whether the angles, mea- 
sured in unit of ir, are rational numbers or irrational numbers. 
For the rational cases, the conductivity diverge as power law 
with exponent /3 = 0.22 and the MSD a = 1.178. While for 
the irrational case, the conductivity saturates to 0.225 and the 
diffusivity to 0.15. 

All of these works in the billiard models lead us to think 
about the intrinsic relation between normal (anomalous) par- 
ticle diffusion and normal (anomalous) heat conduction. Con- 
sequently, there were two independent works trying to relate 
the diffusion exponent a to f3, one by Li and Wang [54] and 
the other by Denisov et al. [75]. Both the studies concen- 
trated on the gas models. Particles transport between the ends 
of dynamic channels and exchange energy with the heat baths 
connected to the ends. In both models, particles do not inter- 
act just as the billiard models described above. The first one 
[54] used a general continuous time random walk description. 
By studying the mean first passage time of the particles, Li 
and Wang finally obtained the relation /3 = 2 — 2 /a. 

In the second study [75], Denisov et al. focused on the spe- 



cial Levy type random walk, which generally assumes that a 
particle waits at each point for a waiting time following dis- 
tribution ip w (t) ~ t~ lw ~ 1 (-f w > 1) and then make a jump 
with distance following distribution ip/(x, i). For Levy walk 
[142], %j) f (t) ~ t-^- 1 andV/OM) ~ ^f(t)S(\x\-vt). This 
means that during a single flight, the particles travel at a con- 
stant speed v. Depending on j t , the particles diffuse according 
to (x 2 (t)) ~ t a with a = 2 for < 7/ < 1; a = 3 - 7/ for 
1 < 7/ < 2 and a = 1 for 7/ > 2. For the superdiffusive 
cases 1 < 7/ < 2, the authors obtained the mean first passage 
time t oc L 7/ . The heat flux contributed from a single particle 
then scales as J oc L~ 7 ^ = L a ~ z . Considering that the tem- 
perature gradient VT oc L^ 1 and the total number of particles 
in the channel N oc L, it was obtained that the heat conduc- 
tivity k oc oc L a ~ x which means j3 = a — 1, Denisov 
et al. also considered Levy flight with jump length distribu- 
tion i/jf(x,t) = ipf(x)5(t — tf) where ipf(x) ~ x~^f . 
However, due to the divergence of mean square displacement, 
(\x\) 2 ~ t a was used instead of the mean square displacement 
which gave j3 = 2/7/ — 1 = 2/(3 — a) — 1. For the sub- 
diffusive cases, the conductivity vanishes because the mean 
passage time diverges. 

Although the billiard gas models helped us clarify a few 
puzzles on the conditions for diffusive heat transport, we 
should admit that due to the lacking of interactions and lo- 
cal thermal equilibrium as pointed out by A. Dhar and D. 
Dhar [143], conclusions drawn from them could not be di- 
rectly transferred to oscillating lattices. The properties of the 
billiard gas models are dominated by their dynamic proper- 
ties, or even simpler, geometric properties, while for lattices, 
thermodynamic properties are more important and of interest. 
Moreover, heat flux in lattices takes place by phonon transport 
without net particle flow, which also make the lattice systems 
different from the billiard gases in underlying physics. 

For energy diffusion in lattices, Cipriani et al. studied a 
one dimensional diatomic hard-point model [144]. The model 
consists of a chain of hard-point particles with alternating 
masses, m<n = m and rri2i+i = rm, lying on a line segment 
with length L. The particles move in one dimension and elas- 
tically collide with other nearest neighbors when they meet at 
the same point. Therefore, the order of the particles are con- 
served and the system is more like a lattice model compared 
to gas models. In their simulation, the average kinetic energy 
is chosen to be (rriivf} /2 = 1 and m = 1, 

It is clear that when r = 1, the collisions only lead to ex- 
change of velocities and the dynamics is integrable. When r 
deviates from 1, the system is non-integrable but still remains 
non-chaotic because the evolution equations is linear: if the 
velocities of two successive particles before collision are Vij, 
then after collision they will change according to a linear form 

v' i = v j ± Y ^ ; (v i -v j ), (14) 

where plus (minus) sign in ± corresponds to i even (odd) and 
j = i±l. 

After introducing an infinitesimal perturbation 5v{ (t = 0) 
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to the i'th particle, the spread of the energy perturbation at 
later time can be characterized by <J( 2 ) (i,t) = rrii(Svi(t)) 2 . 
This quantity can be regarded as the energy distribution in the 
energy diffusion process. It is observed by numerical simu- 
lation that the profile of <W)(£, i) (See Fig. 15) satisfies the 
power-law ansatz 6(2\(i,t) as t r 6( 2 ) (*A 7 ) very well. A best 
fit for i — gives 7 = 0.606 ± 0.008. The evolution of 
the profile <W) (i, t) was then compared with Levy walk with 
waiting time distribution %jj(t) = [142]. The scaling 

law of Levy walk is [ 145] 
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for the central part. Therefore, fj, = 1/7 as 5/3. Numerical 
simulations showed a good agreement between this two mod- 
els except that the ballistic wavefronts are slightly broadened 
(Fig. 15). Since for Levy walk, the mean square displacement 
(a: 2 (i)) oc t a with a = 3 — fi, therefore a should be 4/3 
approximately. Numerical calculation of mean square dis- 
placement of their diatomic gas model showed that f3 as 1.35, 
which supported the connection of their model with Levy 
walk. 

The heat conduction of this model had been separately 
studied earlier by Hatano [5], Dhar [6] and by Grassberger 
et al. [7]. The first work [5] used r = 1.22 and found that 
p = 0.33 ~ 0.37. The second one [6] observed a small diver- 
gence of the heat conductivity with (5 as 0.17 when r = 1.22 
using length N up to 1281. The last one [7] used a very effi- 
cient event driven algorithm and simulated chains with length 
up to 16383. It was found that for very large N, the divergent 
exponent tends to (3 = 0.32^° 01 f° r a ^ mass ratio between 1 
and ~ 5, which was claimed to be consistent with the predic- 
tion of renormalization group theory, /3 = 1/3, by Narayan 
and Ramaswamy [45]. 

The hard-point gas model cannot mimic the real lattice 
models because the interactions between particles are present 
only when they collide. For interacting lattices connected by 
springs, another important systematic study on energy dif- 
fusion were carried out by Zhao [146]. In this work, Zhao 
showed a heuristic approach to investigate the energy distri- 
bution in a diffusion process, denoted as p(i, t), using equilib- 
rium energy spatiotempral correlations, i.e. 
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This quantity is in equivalence with 5(p)(i,i) in Cipriani 
et aZ.'s work [144]. Using this method, Zhao studied three 
well known lattices showing different behaviors of heat con- 
duction: the Toda lattice which is ballistic, the FPU-/3 lat- 
tice which is superdiffusive and the </> 4 lattice which is dif- 
fusive. It was found that for all these models, the energy 
correlation functions can well capture the characteristics of 
the corresponding energy diffusion (Fig. 16). Specifically, 
for Toda lattice, the energy distribution is dominated by two 
clear ballistic wave front; for 4 lattice, the distribution is 




FIG. 15. The rescaled perturbation profiles 5(2) for the HPG 
model at t = 40, 80, 160, 320, 640, 1280, 2560, and 3840 for 
7 = 3/5. Inset: the perturbation profile of at t = 640 (solid 
line) is compared with the propagators of a Levy walk with expo- 
nent p — 5/3 and velocity it = 1 (dotted line). For further details 
seeRef. [144]. 
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FIG. 16. p{i,t) versus t for (a) the Toda lattice, (b) the FPU-/? 
lattice and (c) the </> 4 lattice at time t = 100 (blue), 200 (green) and 
300 (red), respectively. For further details see Ref. [146]. 



Gaussian-like; for FPU-/3 lattice, it is a combination of the 
two above, two seemingly ballistic peak and a Gaussian-like 
central part. Quantitative results for the mean square displace- 
ment (a; 2 (t)) = were calculated, which shows 
power law dependence on time (x 2 (t)^ ~ t a . The exponent 
a as 2 for Toda lattice, a as 1.4 for FPU-/? lattice and a as 1 
for (j> 4 lattice, which coincide with their ballistic, superdiffu- 
sive and diffusive behaviors, respectively. Moreover, it was 
argued that this supports the /3 = a — 1 relation because 
the exponent j3 for FPU-/? model is 2/5 as predicted by the 
Peierls-Boltzmann theory. 

In a recent work, Zaburdaev et al. restudied the energy 
diffusion in one dimensional systems using the phenomeno- 
logical Levy walk approach [147]. They first generalized 
the Levy walk model to involve random fluctuation in the 
movement of particles. Specifically, compared to the standard 
Levy walk, it was additionally assumed that during a single 
flight between two successive collisions, a particle's position 
evolves according to Langevin equation x(t) = vq + £(t), 
where is a Gaussian random process with correlation 
(£(i)£( s )) = D v 6(t — s). When the randomness is absent 
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FIG. 17. Left: rescaled perturbation profiles of the HPG model at 
time t = 1000,2000,4000,6600 respectively. The profiles are scaled 
under Eq. (15). The inset shows the ballistic humps after the scaling 
transformation Eq. (17). Right: energy correlation functions e(i,t) 
for the FPU-/? model (thin solid lines) and the propagators of the 
generalized Levy walk model (thick dashed lines). For further details 
see Ref. [147]. 



(D v = 0), this model reduces to the standard Levy walk 
model with moving speed vq. Analytic analysis showed that 
for their generalized Levy walk model, the central part of the 
particle density profile P(x, t) follows the scaling of standard 
Levy walk, i.e., Eq. (15), while the two ballistic humps follow 
another scaling 
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where x = x — Vot. 

Using this model, Zaburdaev et al. restudied the energy 
diffusion in hard-point particles (HPG) model [5-7, 144] and 
the FPU-,3 model using Eq. (16) by Zhao [146]. For the 
hard-point particles model, Zaburdaev et al. calculated the 
evolution of an infinitesimal perturbation p(i,t) following 
Ref. [144]. It was found that the scaling relations Eq. (15) 
with 7 = 5/3 and Eq. (17) are fulfilled separately for the 
central part and the ballistic part, respectively (Fig. 17). For 
the FPU-/? model, the authors performed a direct compari- 
son between the energy spatiotempral correlation functions 
e(i,t) as in Eq. (16) and the probability density for the gen- 
eralized Levy walk model. The comparison was made at time 
t = 1000 and 2000 respectively (Fig. 17). A good agree- 
ment was observed by choosing the average energy per parti- 
cle e =l,v = 1.384 and D v = 0.49. 

Although there are a few numerical studies showing that 
energy diffusion in some one dimensional systems can be de- 
scribed by (generalized) Levy walk, it can be argued that a 
fundamental theory connecting these two phenomena does not 
exist in general. All existing works are based on phenomeno- 
logical approaches. It is not known either for which kind of 
systems, the Levy walk description is applicable. 

However, regarding the connection between (anomalous) 
heat conduction and (anomalous) energy diffusion, in a most 
recent work, Liu et al. have rigorously proved an equality 



that can be generally applied to all homogeneous one dimen- 
sional systems near thermal equilibrium [148]. In their work, 
Liu et al. considered a general one dimensional Hamiltonian 
system which was initially (t = — oo) at thermal equilibrium 
with temperature T. Then by assuming local thermal equi- 
librium, a small perturbation with perturbative Hamiltonian 
Hext = — h J e(x,t)AT(x)dx was imposed on the system 
from t = — oo, where e(x, t) is the local energy density. This 
form of perturbation will drive the system away from equilib- 
rium and finally reach a steady state with temperature profile 
T(x) — T + AT(x), provided that local thermal equilibrium 
is satisfied. After removing the perturbation at t = 0, the sys- 
tem will start to relax to equilibrium at temperature T again. 
This relaxation process is interpreted as energy diffusion. In 
their work, Liu et al. first analytically proved that in the diffu- 
sion process (t > 0), the normalized distribution of the energy 
profile can be written as 
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where C ££ = (Ae(x, i)Ae(0, 0)) is the autocorrelation func- 
tion of energy density fluctuation. Consequently, it was 
proved that the mean square displacement of energy diffusion 
should satisfy a very general equality independent of the sys- 
tem Hamiltonian and the initial temperature profile AT(x), 
which reads 



d 2 x 2 {t) C(t) 



dt 2 



kBT 2 cv ' 



(19) 



where C(t) is again the heat flux autocorrelation function ap- 
pearing in the Green-Kubo formula Eq. (8). Because in their 
derivation there is no any explicit assumption except the local 
thermal equilibrium, this equality should be universal, even 
for isotropic systems in higher dimensions. 

The authors also showed a few corollaries of this equal- 
ity for different diffusive behaviors. When the system is nor- 
mal diffusive, i.e. x 2 (t) ~ 2Dt in large t limit, by inserting 
Eq. (19) into the Green-Kubo formula, we will get k = cyD. 
For the superdiffusive case x 2 (t) ~ t a with 1 < a < 2, the 
Green-Kubo formula with a cut-off time t c oc L will lead to 
k ~ L"" 1 , which means /? = a — 1. While for subdiffu- 
sive cases (0 < a < 1), the heat conductivity decays to as 
k - l/L 1 -" [148]. 



CONCLUSIONS 

In this Colloquium we took the reader on a tour present- 
ing the state of the art of the topic about anomalous thermal 
conduction and heat diffusion in low dimensional nanoscale 
structures, an emerging research direction which is expected 
to shade light for renewable energy and thermal management 
for future nano devices. Particularly, we surveyed and ex- 
plained physical mechanisms that respond to the anomalous 
thermal conductivity. 
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We firstly reviewed recent experiments on thermal conduc- 
tivity of nanotube, nanowire and graphene. Then, we pro- 
vided a review on the computational investigations of ther- 
mal conduction in nanascale, concentrating on atomistic sim- 
ulations. At last, various theoretical models to explain the 
anomalous thermal conductivity and the connection between 
thermal conduction and heat diffusion were overviewed. Both 
numerical and experimental studies in nanostructures rang- 
ing from nanotubes [8, 10, 51] nanowires [9], to polyethy- 
lene nanofibers [17, 26, 149] have shown that heat con- 
duction in low-dimensional systems does not obey Fourier's 
law even though the system length is much longer than the 
phonon mean free path. In these structures phonons trans- 
ports superdiffusively — a process faster than random walk but 
slower than ballistic motion, which leads to a length depen- 
dent thermal conductivity. Thus nano structured materials are 
promising platforms to testify fundamental phonon transport 
theories. 

Moreover, given the fact that thermal management de- 
vices generally can be considered for applications that require 
power ranging from milli-watts up to several thousand watts, 
the general range of applications where they are indispensable 
is stupendous. Compared with decades ago, we have a better 
and clearer understanding of the thermal conduction and heat 
diffusion in nanoscale materials. 



the interface by taking into account the anomalous thermal 
transport characteristics of nanostructures. To provide quan- 
titative prediction from theory and simulation, how to set up 
a transport theory by incorporating nonlinearity in the quan- 
tum regime is still a challenge. The nonequilibrium Green's 
(NEGF) function method [150] certainly serves as an ele- 
gant mathematical framework. However, when including the 
phonon-phonon interactions, the NEGF presents a cumber- 
some task. Thus further systematical investigations combine 
experimental and theoretical efforts on fundamental mecha- 
nisms will be helpful to advance the field. 
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OUTLOOK AND CHALLENGES 

However, there are still many open questions. On the exper- 
imental side, experimental demonstration of low-dimensional 
phonon transport theory depends on how to measure its ther- 
mal conductivity accurately. To this end, one needs to make 
the sample smaller, for example towards a dozens of nanome- 
ters or even a few nanometers only. Difficulties in integrating 
nanomaterials with suspended structures suitable for thermal 
measurements have resulted in rare experimental studies of 
size-dependent thermal conductivity in ID and 2D systems. 
And another challenge is to get rid of the contact thermal re- 
sistance. The unknown thermal contact resistance remains the 
primary technical challenge in testing the Fourier's law. More 
efforts are needed to explore a brand new method of measur- 
ing thermal contact resistance directly. 

Because of the complexities inherent in thermal contact re- 
sistance measurement, this task needs the help from theoreti- 
cal predictions. However, so far a truly comprehensive theory 
for this effect is lacking. The current approaches for ther- 
mal transport across an interface, such as the acoustic mis- 
match (AMM) theory and the diffusive mismatch (DMM) the- 
ory, are based on the assumption that phonon transport via 
a combination of either ballistic or diffusive transport on ei- 
ther side of the interface. Both schemes offer limited accu- 
racy for nanoscale interfacial resistance predictions, due to as 
we have shown above, phonons transport super-diffusively in 
low-dimensional systems. Therefore, it is necessary to estab- 
lish an improved theory describing thermal transport across 
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